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Abstract 

The symbol of the Hessian for an aeroelastic optimization model problem is ana- 
lyzed. The flow is modeled by the small-disturbance full potential equation and the 
structure is modeled by an isotropic (von Karman) plate equation. The cost function 
consists of both aerodynamic and structural terms. In the new analysis the symbol of 
the cost function Hessian near the minimum is computed. The result indicates that 
under some conditions, which are likely fulfilled in most applications, the system is 
decoupled for the non-smooth components. The result also shows that the structure 
part in the Hessian is well-conditioned while the aerodynamic part is ill-conditioned. 
Applications of the result to optimization strategies are discussed. 
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1 Introduction 


Lately, there is a growing interest in Multidisciplinary Design and Optimization (MDO) 
[l]-[4]. An important problem in that field is the aeroelastic optimal design problem (for 
example, [5]- [7]). In this problem there are two coupled disciplines: aerodynamics and 
structural analysis. The problem is to compute the aerodynamic shape and structural rigidity 
such that some given cost function is minimized. 

The purpose of this paper is to demonstrate new analysis of Hessians for MDO problems 
on the above aeroelastic optimization problem and to draw some practical conclusions. The 
approach is to consider a simple model problem and compute the symbol of the Hessian near 
the minimum for the non-smooth frequencies. The Hessian contains curvature information 
which is essential for the solution of ill-conditioned optimization problems. Hessian symbols 
were previously computed for smoothing predictions in the development of multigrid one-shot 
methods [8]- [1 1] and lately for the analysis of inviscid aerodynamic optimization problems 
[12]. The analysis in this paper indicates that for the non-smooth components the system 
is decoupled under certain conditions, which are likely fulfilled in most applications. The 
analysis also shows that the structures part in the Hessian is well-conditioned while the 
aerodynamics part is ill-conditioned. 

One consequence of this result is that if the decoupling condition holds the solution of 
such problems can be achieved in two stages. In the first stage, the MDO approach should 
be taken on a coarse model; that is, the flow and the structure equations are considered 
simultaneously during the minimization, which is a more complex problem than optimizing 
the decoupled individual disciplines problems. In the second stage, a refined CFD code for 
the flow and a detailed finite element code for structure should be used in a serial algorithm in 
which the shape is optimized relative to aerodynamic considerations, followed by structural 
optimization limited to a given shape that under flow conditions it will twist and bend to the 
aerodynamic optimal shape [13, 14]. This approach should result in a good approximation 
of the multidisciplinary optimal solution. 

A second consequence is that if the decoupling condition does not hold, and thus the two 
disciplines are coupled, a multidisciplinary algorithm is necessary and a preconditioner for 
the coupled system is necessary to obtain effective convergence. 

The paper outline is as follows. In Section 2 the optimization problem is formulated. 
In Section 3 the necessary conditions for a minimum are derived with the adjoint method 
and their relation with the Hessian is discussed. In Section 4 the symbol of the Hessian for 
the non-smooth frequencies is derived by using local mode analysis. In Section 5 decoupling 
conditions and multidisciplinary preconditioners are derived. In Section 6 applications of 
the result are finally discussed. 


2 Problem Formulation 

In this section the aeroelastic analysis problem and the optimal design problem are presented. 
The aeroelastic analysis problem couples the full potential flow equation with the isotropic 
von Karman plate equation to give the pressure distribution over the plate, p, and the 
plate deformation, W, for a given plate shape, a, and rigidity distribution, D. The design 
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problem is to compute the “best” shape and structural rigidity so that a given cost function 
is minimized. 

The cost function is composed of aerodynamic and structure parts. The aerodynamic 
cost function estimates performance by measuring the difference, in L 2 norm, of the pressure 
distribution from a desired one. The structure cost function gives a measure of the structural 
weight and penalizes structural deformation. 

Since our interest is in a local mode analysis of the Hessian near the minimum, we consider 
the small disturbance equations of flow over a flat plate. 

2.1 The Flow Model 

We choose the full potential equation as a model for the flow. It approximates inviscid 
flow characteristics and is used in applications for aerodynamic optimal design (for example, 
[15]). For the analysis of the cost function’s Hessian in the vicinity of the minimum it is 
enough to consider small perturbations of the shape from the optimal solution. The resulting 
changes in the potential satisfy the steady state small disturbance full potential equation. 
The geometry is taken to be half-space 0 = (x,y,z > 0), where the x axis is the stream- 
wise coordinate, y is the coordinate perpendicular to the stream and parallel to the plate 
(spanwise direction), and z is in the normal direction to the plate. 

The Aerodynamic State Equation: 


L <f> = 0 z > 0 

B <j> = C(a + W) z = 0 

with the following definitions of the interior operator, L, and boundary operators, B and C, 

L = (1 — M^)d X x + d yy + d zz 
b = a 
c = Uoodx 

where {/<*, is the free stream velocity, M is the Mach number, with the following far-field 
boundary conditions: 

Inflow Boundary Condition 
Subsonic: <f> x = Uoo 

Supersonic: <p x = Uoo and <f> = 0 (we assume that the normal free stream velocity, 14o, is 
zero). 

Outflow Boundary Condition 
Subsonic: <f> x = U^ 

Supersonic: No Boundary Condition. 

The missing low-order terms in the boundary condition of (2.1) vanish if the analysis is 
performed around a flat shape. 


( 2 . 2 ) 
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2.2 The Structural Model 

The structural model consists of the isotropic von Karman plate equation for the displace- 
ment W [16, 17]: 

G{D,W) = pC<t> z = 0 (2.3) 

with the following definition of the the operator G: 

G(£>, W) = ( DW XX ) XX + ( DW yy ) yy + v[{DW yy ) xx + {DW xx ) yy \ + 2(1 - v){DW xy ) xy , (2.4) 

where D is the plate rigidity distribution, p is the flow density and v is the Poisson ratio. In 
two space dimensions Eq.(2.3) reduces to the beam equation: 

( DW XX ) XX = pUootx z = 0. (2.5) 

The boundary conditions for the plate are that it is clamped at y = 0 (or at x = 0 for a 
beam) and free at the other boundaries: 

W(x, 0) = W y {x, 0) = 0. 

However, Eq.(2.3) is elliptic, so the effect of a high-frequency error component in the deflec- 
tion W is local, and therefore the plate boundary conditions do not play a role in the local 
mode analysis. 


2.3 The Cost Function Model 

The definition of the cost function is not unique and depends on the specific application 
under consideration. In general the requirement of the aeroelastic optimal design is that it 
have maximum aerodynamic performance and minimum structural weight and deformation. 
Some of the desired features of the final design are in many cases modeled by a set of 
inequality constraints, as is the case for the minimum deformation requirement. However, 
for the purpose of this paper we will avoid inequality constraints by adding a term to the 
cost function which penalizes deformation. In the following the different terms composing 
the cost function are discussed. 


• The Aerodynamic Performance Term 

A common aerodynamic cost function is drag (or drag over lift). However, in inviscid 
aerodynamic optimization models a commonly used cost function is pressure matching (for 
example, [18]- [22]). In potential models the pressure, p, is related to the potential, <f>, by the 
Bernoulli relation 

P — pU( oo^x- 


This results in the following cost function term 

F aeTO = - rfda 

where da is an integration element on the shape T. The target distribution, f*(x, y ) 6 L 2 (T), 
is related to the desired pressure distribution, p*(x,y), by the relation 




fj^y) 

pUoo 
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• The Structural Weight Term 

Another important factor in aeroelastic design is the resulting weight of the structure. In 
practice the weight is measured by the sum of the weights of all the components composing 
the structure. In plate models the weight is related with the plate rigidity, D , and is given in 
the following table, where E is the Young modulus of elasticity, b and h are the cross section 



D 

weight 

beam 

Ebh J 

p p f r bhdx 

plate 

eP 

12(1— i/ 2 ) 

p v f r tdxdy 


components of the beam, p p is the structural density, and t is the plate’s thickness. 

In both cases the weight of the structure is proportional to where d is the space 
dimension: 

pweight ^ j D \ d(J 

• The Structural Deformation Term 

As a result of the pressure, p, exerted on the plate by the flow, the structure will deform 
its shape by W (bend and twist). In practice the structure is designed so that the amount of 
deformation will be constrained not to exceed some given limits. In this model we account 
for this requirement by penalizing the deformation which is measured by the work of the 
aerodynamic pressure on the plate, pW . This will add to the cost function the following 
term (using the Bernoulli relation): 

F deS ° rm = pU^ J^ X W do. 

2.4 The Optimization Problem 

We define the cost function, F = F(<j > , W, D ), to be 

F(<j>, W, D) = J^{<j> x - f*) 2 do + 72 ^ D*do + 73 J^<f> x Wdo (2.6) 

where 71,72,73 are parameters. The cost function is a map from a function space to IR. 

The minimization problem is to find a shape function, a, and rigidity distribution, D , 
such that the cost function is minimized subject to Eqs.(2.1) and ( 2 . 3 ). We assume the 
existence of a solution for both the state equations and for the optimization problem (a 
rigorous treatment of existence and uniqueness of solutions is beyond the scope of this 
paper) . 


3 Adjoint Formulation and the Hessian 

In this section the necessary conditions for a minimum are derived with the adjoint method 
[l 9 ]-[ 24 ]. The necessary conditions are given as a set of state equations (the analysis prob- 
lem), costate equations (the adjoint problem) and design equations (optimality conditions). 
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Then the relation between the design equation residuals and the Hessian of the cost function 
is discussed. This relation will be used in the next section to derive the Hessian’s symbol. 


3.1 The Necessary Conditions for a Minimum 

The Lagrangian is a functional defined by 

C(<f>, W, a, D, £, A, rj) = F(<f>, W, D) + f r £ (B <j> - C (a + W))da+ 
f Q XL<f>dQ + f r t] (G(D, W) - pC(f>)da 

where £ = rj = rj(x,y) and A = \(x,y,z) are Lagrange multipliers. The first order 

necessary conditions for a minimum are derived by the requirement that the first order 
variation of the Lagrangian vanish (this is known as the adjoint method and the resulting 
conditions are known as the Kuhn-Tucker conditions). 

When considering the variation of the structure state equation a linearization is per- 
formed, 

G (D* + D, W* + W) = G (£>*, W m ) + G b (W*)D + G w {D*)W + h.o.t. (3.2) 

where D and W are small perturbations of the displacement and rigidity from the optimal 
solution W* and D * respectively, and where the linearized operators Gd and Gyv are defined 
as follows 

g d(w-)d = D xx w; x + D yy w; y + */ [b xx w; y + b yy w; x ] + 2(1 - v)b xy w; y 

Gw(D*)W = G{D\W). 

Formally, W* and D * serve as non-constant coefficients in the linearized structure operator. 
In the following the costate and design equations are given. 


Costate Equations 

LA = 0 z>0 

BA + pCy = F 4 , z = 0 
G\v(D*)rj + CA = — Fw 2 = 0 

Inflow Boundary Condition 
Subsonic: A* = 0 

Supersonic: No Boundary Condition 

Outflow Boundary Condition 
Subsonic: \ x = 0 
Supersonic: A = 0 and A x = 0 


Design Equations 


CA = 0 2 = 0 

G d (VF*)t? + F d = 0 2 = 0 


(3.4) 


(3.5) 
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where 


F<t> = -271 { 4 > x - f*)x - izW x 

Fw = 73<^x (3.6) 

F d = *D^, 

and where the operators in the adjoint and design equations (3.4-3. 5) satisfy 


L = L 

C = -C 
G w (D*) = G W (J D*) 
G D (W*) = Gd(VT*). 


The adjoint boundary operator B corresponds to the normal derivative, d z , applied to a 
solution of the interior costate PDE, A, when using the adjoint far-field boundary conditions. 
We assume the existence of a solution to the costate equations. 


3.2 The Relation of the Hessian with the Necessary Conditions 

If the state and costate equations are satisfied, in the strong sense, then the variation of the 
Lagrangian (3.1) is equal to the variation in the cost function and is given by 

SC = J aC\d(T + j D(G J >(W m )ri + F D )d(j (3.7) 

where o; and D are variations in the design variables. Therefore the quantities multiplying 
a and D in (3.7) are the sensitivity gradients of the cost function with respect to the design 
variables, when computed on the constraint manifold: 

Ji = CA 

52 = Gd(W*) 77 + Fp. 


(3.8) 


The state and costate equations, (2.1), (2.3) and (3.4), give an implicit relation between the 
costate variables and the design variables: 

A = A(or, D) 

v = 

Using equations (3.8) and (3.9) we can write the following relation which holds near the 
minimum (a* and D* denote the optimal value of the design variables a and D, respectively): 

<7i (A (a* + a, D* + D)) = H\\6t + H^D + h.o.t. , . 

g 2 {D~ + D, via* + a,D* + D )) = H 21 a + H 22 D + h.o.t. [ } 

where at the minimum 


J1 <A(o-,D , ))= SJ (fl%,(a*,fl*)) = 0. 

We conclude that on the constraint manifold, near the minimum, the Hessian of the cost 
function relates the errors in the design variables with the residuals of the design equations 
(sensitivity gradients). In the next section we will use this fact to calculate the symbol of 
the Hessian. 
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4 Derivation of the Hessian’s Symbol 

In the following section we compute the symbol of the Hessian with local mode analysis. 
Hessian symbols were previously computed for smoothing prediction in the development 
of multigrid one-shot method [8]- [11] and lately for the analysis of inviscid aerodynamic 
optimization problems [12]. In the following the local mode analysis is outlined. 

The analysis is performed in the vicinity of the minimum where the design variables are 
assumed to have an error ct and D. We assume that the state and costate equations are 
satisfied and consider the errors in the state and costate variables ( 4>,W,\,fj ) with respect 
to their value at the optimal solution. These errors are assumed to satisfy homogeneous 
equations similar to Eqs.(2. 1,3.4, 3.5), and a linearization of Eq.(2.3). We then consider the 
high-frequency errors in the design variables and compute an explicit solution of the problem 
in terms of exponential functions in a half-space. Then with a standard procedure the 
problem in a half-space is reduced to the boundary. On the boundary we study the mapping 
from the transformed design variables errors to the residuals of the design equations, (fiq, (fe)- 
The symbol of this mapping gives the eigenvalues of the Hessian. 


4.1 Fourier Representation 

We start with the Fourier representation of the solution in a half-space and perform local 
mode analysis. Consider errors in the solution of the form 

S(x,y) = d(w 1> a> J )e i <"^»> , 41 , 

D(x,y) = 0(u 1 ,u,)e i <'-' 4 “">. 

As a result, the errors in the state and costate variables are assumed to have the following 
form: 

t(x,y,z) = 

A (*,»,*) = A( Wl ,cu 2 ,cu 3 )e‘^^^3d 

W{x,y) = W(ui, 
fj{x,y) = fjy 

Before computing the relation between the state and costate error symbols, (</>, A, W,rj), and 
the design error symbols, (a,D), we reduce the problem to the boundary by eliminating u> 3 
from the symbols <p and A. 

4.2 Reduction to the Boundary 

The reduction to the boundary is done by eliminating u> 3 from the symbol expressions using 
the interior equations. The following discussion regarding the choice of u> 3 was done in [12] 
and is repeated here for completeness. 

The term <f> satisfies the interior equation for (j>: 

, cu 2 , cu 3 )e'^ lX+u ’ 2y+W3z) =0. (4.3) 
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Assuming a non-trivial solution, 7 ^ 0, Eq.(4.3) results in an algebraic relation between u> x , 
u >2 and u> 3 : 

(1 - M 2 )u\ + <J 2 + (4 = 0. (4.4) 

The choice of u> 3 should be done such that it will result in a physical solution. We differentiate 
between subsonic and supersonic flows. 

4.2.1 Subsonic Flow 

In the subsonic regime (M < 1) the physical solution is given by 

u 3 = iyjw\( 1 - M 2 ) +u>2, 

which corresponds to decaying solutions: 

4>(x,y,z) = 4>(u x , W2 )e i ( w ‘* + "»»)c-(V'^( 1 - M3 )+^)* 

\{x,y,z) = A( Wl , u 2 )e' {uJ ' x+w * y) 

In that case the symbols of the boundary operators, B and B, are given by (recall that B 
and B are the normal derivatives applied to a solution of the interior state and costate PDE 
respectively) 

B = B = -^(1 -M 2 )+u;|. (4.5) 


4.2.2 Supersonic Flow 

We differentiate between two supersonic cases which are determined by a flow speed denoted 
M c and given by 



The case 1 < M < M c results in the same symbols for B and B as for the subsonic flow case 
(Eq.(4.5)). 

In the case M c < M both signs of u 3 in (4.4) correspond to physical solutions. The 
positive root correspond to the characteristic which propagates into the shape, <f>+, and the 
negative root correspond to the characteristic which propagates out of the shape, <£>_, (and 
a similar expression for A): 

^(s,y,z) = ^> + {u Xl u>2) X+W2 y+ V / I w i ( 1 --W 2 )+“! I 2 ) + <j,_( UuU2 ) c *'(" i*+w- VH l 1 W . 

( 4 - 6 ) 

Since the inflow information does not change as a result of a shape perturbation, the 
following holds: 

4>+(x,y,z) = 0. (4.7) 

In the same manner the outflow characteristic of the adjoint is not changing as a result 
of a shape perturbation: 

A_(i,y,z) = 0. (4.8) 
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Therefore 


4>(x,y,z) = 4 >-(u u lj. 2 ) *+^y~V K C 1 - ^ 2 )+- 2 2 h) 

A (x,y,z) = A + (^ 1 ,u; 2 )e t(l " ir+ ‘" 2j ' + ^ (1 - M2)+ ^ |2) . 

We conclude that for flow speeds M c < M the boundary operator B is antisymmetric, 
(with respect to the adjoint operation), and the symbols B and B are given by 

B = - 

B = i 

In all flow conditions the multiplication BB results in the same expression: 

BB = |o^ (1 - M 2 ) + w||. (4.10) 

By eliminating u >3 from the transformed equations the state and costate flow equations can 
be written on the surface (cjj ,uj 2 ) which corresponds to the boundary (x, y). 

4.3 Treatment of the Structure Equations 

In this subsection we give a short note concerning the transformation of the structure state 
and costate equations. The structure state and costate equations contain non-constant 
coefficients which should be frozen prior to the local mode analysis. The structure state and 
costate error equations are given by (see Eqs.(2.3,3.2,3.4)) 

G d (W*)D + G w (D*)W = pC4> z = 0 f4n) 

Gw(f}*)?j T C A = —Fw z = 0 

where Z), W, 4>, ij and A denote the error variables, Fw denotes the error in Fw, and the 
operators Gjj and G\v are defined in (3.3). Since Eqs.(4.11) have variable coefficients, D 
and W*, it is necessary to freeze them around a point on the boundary. This procedure is 
justified as long as the errors in the design variables are highly oscillatory compared to W * 
and D*. We denote the values of W*(x,y) and D*(x,y) at a point (x 0 ,y 0 ) on the boundary 
by Wq and D respectively. 

4.4 The Coupled State and Costate Equations in Fourier Space 

In terms of their Fourier representation on the boundary, the state and costate error equations 
are given by the following matrix form: 

( B -C 0 0 \ ( 4> \ ( Ca 

-pc G w {Do) 0 0 w -G d (W*)D 

—F^q —F^w B pC A 0 

y Fw<)> 0 C Gw{D'o) J \ V J V 0 
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The various symbols are given explicitly by 


= 2 7i 

F^w = — *73 w i 
Fw4> = ^73^1 

Gw(T>S) = I> 5K 2 + u ; 2 2 ) 2 + /.oJ. 

Gd(w 0 *) = + .H^) - u 4 (W^, + vw;„) - ^(2(1 - v)W^ y ) 

C — ^ u M . 


(4.13) 


Note that the terms originating in the cost function serve as a coupling symmetric block 
between the state and costate systems. 


4.5 The Symbol of the Hessian 

The design equations residuals, in the transformed space, are given by 


h = C\ 

92 = d D (WS)rj + F dd (D- 0 )D 


(4.14) 


where 


Fdd 


72(1 ~ d) 

d 2 



3 —2d 
d 


and the symbols <71 and 92 are the symbols of the sensitivity gradients <71 and <72 , respectively 
(see (3.8)). 

Inverting the system (4.12) and substituting A and 17 in the symbol of the design residuals 
(4.14) results in a relation between the residuals of the design equations and the errors in 
the design variables. In Fourier space, 


( 9i \ ( Hn 

\h ) \ H n 



(4.15) 


where the matrix H X3 is the symbol of the Hessian, as discussed in Sec. 3.2. Hu is the symbol 
of the aerodynamic optimization Hessian, H 22 of the structural optimization Hessian, and 
#12, H 21 are the coupling terms. In the following, the terms H{j are given explicitly: 


£ = -C 2 Gw(F^Gw + 2pCFj W ) 

( BG W - C 2 P ){BGw - C 2 P ) 

f. _ CGp{CP^Gw ± BF^wGw ± pC^w) 
( BG W - C 2 p)(BG w - C 2 p) 

£j _ CG d {CF u G w + BF^wGw + pC 2 F^w) 
( BG W - C 2 p)(BG w - C 2 p) 

A CG 2 D {F U C + {B + B)F, W ) , . 

“22 — — : — : 1 - r DE>. 

( BG W - C 2 p){BG w - C 2 p) 


(4.16) 

(4.17) 

(4.18) 

(4.19) 
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Since Gw is a fourth order polynomial in u>i, 2 , Gp, F u , B and B are of second order, F^w 
and C are first order, and Fd is of zero order, the principal parts of the above expressions 
(the asymptotic limits of high-frequencies), are given by 


H 


-C 2 F, 


4><t> 


11 


BB 


— 271 ul 


u. 


|w, 2 (l -M 2 )+u> 2 2 I 


(4.20) 


H l2 = H n 


C 2 GdF$<p 

BBG W 


_ + ■'Wo-,,) + U^(W^, + vWjLj + 2a,, ^(1 - v)W^) 

Dl |w?(l - M 2 ) + • (uj 2 + u 2 ) 2 

(4.21) 


Hn » F dd = ^ (D’ 0 ) ““ . (4.22) 

Note that for simplicity we assumed a complex representation of the errors, (4.1), and 
obtained a complex Hessian symbol. If considering a real representation, i.e., 

a(x,y) = a(w,, u 2 )e i( -“' x +“* y) + d conj (u; a ,a; 2 )e-‘'( wl ^ 3 ' ) 

where a conj is the complex conjugate of a. and a similar expression for D, then the resulting 
Hessian symbol is real and symmetric. 


4.6 Discretization and the Condition Number 


In practice the problem is solved numerically and thus discretization is introduced. Therefore 
the analysis should be performed in the discrete space and the Hessian will depend on the 
specific discretization. For the “ideal” discretization, the symbol of the Hessian is equal to 
the differential one with the substitution 


(«i,wa) = ( 


0i 

h / 



5 


where (h u h 2 ) are the mesh sizes in the (x,y) directions, respectively, and where and 0 2 
vary in the domain [ — tt, 7t). 

Note that “high-frequencies” are those which obey > c for some constant, c, which is 
determined by the different parameters in the problem. In the discrete space this corresponds 
to 0 t > chi. Since the constant c is independent of the mesh-size h, as the grid is refined the 
portion of high-frequencies in the spectrum increases and therefore the approximation taken 
by the local mode analysis above is more accurate for a larger part of the spectrum. This is 
not surprising since as the grid is refined its resolution increases while the resolution of the 
smooth components remains unchanged. 

The maximum eigenvalue of each of the disciplinary Hessians is estimated by 




Unfortunately the lowest eigenvalue cannot be estimated by the procedure above since this 
is precisely the spectrum range in which the approximation taken by the local mode analysis 
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does not hold. However, the lowest eigenvalue is a fixed number, independent of the mesh- 
size, and therefore the condition number of the Hessian is proportional to A mor . For a 
two-dimensional flow over a beam, (u> 2 = 0), we get for the aerodynamic part of the Hessian 
(see Eq.(4.20)) 

, 271 VI** 1 

|1 — M 2 \ h 2 

We conclude that the aerodynamic part of the Hessian is ill-conditioned and its condition 
number is increasing as the grid is refined (see [12] for further discussion). The structure’s 
symbol (4.22) approaches a constant, for the high-frequencies, independent of the mesh-size. 
We therefore conclude that the structural optimization problem is well-conditioned. 


5 On the Coupling Between Aerodynamic and Struc- 
tural Design 

In the previous section we computed explicitly the Hessian’s symbol. The coupling between 
the two disciplines, during optimization, is determined both by the off-diagonal terms in 
the symbol of the Hessian and by the smoothness of the design variables. We say that the 
optimization problem can be decoupled if (see Eq.(4.15)): 

|£na| > \H 12 D\ (5.1) 

(or equivalently if \H 22 g\\ \H\ 2 g 2 \). If condition (5.1) holds then the aerodynamic opti- 

mization problem can be solved at the first stage, independent of the structural optimization 
problem, setting the structural deformation to zero ( W = 0). Then, at the second stage, 
the resulting shape and potential are given as inputs to the structural optimization problem 
(see more details in the discussion). The solution of this two-stage approach will be a good 
approximation of the multidisciplinary solution, in the high-frequencies. We emphasize that 
if condition (5.1) does not hold then this procedure will result in a poor approximation of 
the multidisciplinary optimal solution. In that case a multidisciplinary preconditioner 
should be constructed to give a “corrected” disciplinary search directions. This is done by 
transforming Eq.(4.15) back to the PDE level. Note that if \H 2 \6t\ \H 22 D\ holds and 

condition (5.1) does not hold, the process of first computing an optimal structure and then 
computing an optimal shape is not well defined. In the following we examine condition (5.1) 
explicitly for two- and three-dimensional flows and illustrate the derivation of a precondi- 
tioner for a simple case. 


5.1 Two Space Dimensions 

In a two dimensional flow over a beam the principal part of the Hessian is given by (see 
Eqs.(4.20)-(4.22)) 


H (ui > 1) = 2 


jz 

li-A^pi 

W'„ 

1 1 — a/ 2 | d; 


hEL w q ~ 
|i— m 2 \ d; 


(5.2) 


12 



In practice the curvature of the deflection in the stream-wise direction, Wq xx , is negligible 
with respect to the rigidity, £>„, as is the case for airfoils. Setting Wq xx = 0 in (5.2) results in 
a diagonal matrix which implies a complete decoupling between aerodynamic and structure 
design. 

In case Wq xx is not negligible, let us assume the optimal solution exists and belongs to 
the following spaces 

aeC s (r) and D E C k {T) 
where s and k are integers. The matrix (5.2) implies 

Hncxocu 2 -* ( 5 . 3 ) 

H n D oc u~ k 

and therefore the decoupling condition (5.1) is satisfied if 

k > s — 2 (5-4) 

i.e., as long as the rigidity D is not much more “rough” than the shape a (D E C s 2 (F) or 
rougher). 

In case condition (5.4) is not satisfied a preconditioner should be constructed. The 
disciplinary search directions, a and D, are solutions of a PDE which is obtained by the 
inverse transform of Eq.(4.15): 


—aca xx — b 2 a = cg x - bg 2 
—acD xx - b 2 D - —bg\ - a(g 2 ) xx 


(5.5) 


with 


a = 


|1-M 2 | 


b = 




w; 


|1-M 2 | £>* 


c= -t72 D* 


Using these directions as the search directions, in a multidisciplinary optimization algorithm, 
should result in a much more effective convergence properties. 


5.2 Three Space Dimensions 

In a three dimensional configuration we differentiate between the stream-wise and spanwise 
directions. Let us assume that the curvature of the deflection in the stream-wise direction 
is negligible, i.e. set Wq xx = 0. As a result the coupling term Hu has the following form: 


Hu{W' 0xx = 0 ) 


Do 




Oyy + u ]Woyy 


+ 2WiU>2(1 — v )^Oxy) 


\u 2 (l - M 2 )+CU 2 | 


(w? + W2) 2 


(5.6) 


For the design of the structure in the spanwise direction only, i.e., freezing the stream-wise 
design as done in practice for aircraft wing design, the off-diagonal terms in the Hessian 
vanish, (u>i = 0), and therefore the problem is decoupled. 
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For errors in the stream-wise direction only, (i.e. u; 2 = 0), the off-diagonal terms in the 
Hessian reduce to 


^12 2 — 0 ) ~ 


2-n 

Oo|l-M ! l' 


By a similar argument to the two-dimensional case, if a € C*(r) and D € C' s_ 1 (r) or 
smoother then the problem can be decoupled. If the structure design variables belong to a 
much rougher space than the shape variables, ( D 6 C s ~ 2 or rougher), then a preconditioner 
should be applied similar to Eq.(5.5) with 


Q = WL 

u |1-M 2 
_ 2^Ul 




° ~~ |1— M 2 | D* 

c=-|72^-t. 

The case a 6 C s and D £ C s-2 needs a more careful examination. In this special case, 
the decoupling condition (5.1) implies 


£»| < \Hu 




Do 


< 1. 


(5.7) 


Let us assume a wing like geometry where a plate of length L is clamped at ( y = 0) and free 
at the other boundaries. Assume that at the optimal solution the plate is bent towards the 
tip (y = L) as a quadratic in y, 


W(x,y) = W u ’!L, 


(5.8) 


where W tip is the maximum deflection at the tip of the plate. In that case, the decoupling 
condition (5.7) becomes 


24(1 - v 2 )vW tip 

WT 2 


< 1 , 


(5.9) 


where E is the Young modulus of elasticity, v is the Poisson ratio and t is the plate thick- 
ness. In case condition (5.9) is not satisfied then the problem can not be decoupled and a 
preconditioner is required. 


6 Discussion and Concluding Remarks 

The symbol of the Hessian for aeroelastic optimization model problem was computed (Eqs.(4.16- 
4.19)). The result indicates that for the non-smooth components the system is decoupled 
under certain conditions. In two dimensions it is enough to assume that the curvature of 
the deflection is negligible, as is the case for airfoils, to obtain decoupling (see Sec. 5.1 for 
further discussion of the decoupling condition). In three dimensions this requirement is not 
enough and unless the rigidity is not much rougher than the shape the system is coupled 
(see Sec. 5. 2). 

The result also shows that the aerodynamic optimization problem is ill-conditioned, and 
therefore second order information is essential for efficiently solving this part of the problem 
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[12], while the structural optimization problem is well conditioned. Thus, it is anticipated 
that the number of optimization iterations required to solve the multidisciplinary problem 
is determined by the aerodynamic optimization part of the problem. 

We now discuss the application of this result to optimization strategies for the solution 
of the problem. We differentiate between two basic approaches, the “disciplinary” and the 
“multidisciplinary”. In the disciplinary approach the solution of the problem is divided so 
that one discipline optimization problem is solved at each stage, decoupled from the other 
discipline. In the multidisciplinary approach both the analysis and optimization solutions 
are performed in a tightly coupled manner. These two approaches are now presented in more 
detail. 

• The Disciplinary Approach - Weak Coupling 

A common practical strategy used to solve large aeroelastic shape optimization prob- 
lems is the disciplinary approach, i.e., design the aerodynamic optimal shape to give 
the best performance and then design a minimum weight structure, restricted to the 
aerodynamic shape, that under flow will twist and bend to the aerodynamic optimal 
shape [13, 14]: 

The Disciplinary Algorithm 

1 . The aerodynamic shape optimization problem is solved setting W = 0, 

min a f r (<f>: t - /* ) 2 da 
subject to 

L<p = 0 z>0 

B<p = Co z = 0. 

2. The structure minimum weight problem is solved given a fixed 
shape and pressure (potential) distribution, i.e. 

mine 7 f r D*da + 7 ' f r <f) x Wdo 

subject to 

G(D,W)=pC<t> z = 0 

where 7 and 7 ' are parameters. 


3. The final shape, /? , is computed such that under cruise 
conditions, (given pressure distribution - p) , the shape will 
deform into the aerodynamic optimal shape o. 
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• The Multidisciplinary Approach - Tight Coupling 

Lately there has been an effort to develop new optimization strategies which couple 
the two disciplines tightly during the analysis and optimization computation. This is 
known as the MDO approach [l]-[7]. According to this approach after each call to the 
optimizer the analysis and adjoint equations are relaxed, or solved exactly, depending 
on the feasibility choice (Multi-Disciplinary Feasibility (MDF), Individual Discipline 
Feasibility (IDF) or All at Once (AAO), [3]). 

The MDO Algorithm 


The coupled aerodynamic shape and structure minimum weight optimization 
problems are solved simultaneously. In order to achieve efficient conve- 
rgence both disciplines gradients play a role in each disciplinary 
optimizer (see Eq.(5.5)), 

min 0vD J r (<f> x - /*) 2 dx + 72 J r D*dx + 73 f r <j> x W dx 

subject to 

L<f> = 0 z > 0 

B<t> = C{a+W ) z = 0 

G(D,W) = pC<f> z = 0. 

The aim of the MDO approach is to couple a refined CFD code with a detailed finite- 
element structural analysis code to compute the aeroelastic states prior to each optimization 
iteration. The computational complexity of the MDO algorithm is much greater than that 
of the disciplinary algorithm since at each multidisciplinary iteration both the aerodynamic 
and structural optimization problems have to be solved (the MDO approach also introduces 
a technical difficulty of joining together two large codes). The result of this paper shows that 
under certain conditions, given in Sec. 5, which are most likely met for aircraft wing design, 
the system is decoupled for the non-smooth frequencies. Therefore, under these conditions 
the MDO approach applied on a fine scale model might not be necessary to obtain a good 
approximation of the optimal solution. The effect of the smooth components can be captured 
by a coarse model containing a relatively small number of design variables and thus can be 
solved by the MDO approach with a relatively low computational cost. This will require 
simple models for the flow (panel method or small disturbance full potential on a coarse 
grid) coupled with a plate model, or coarse finite-element model, for the structure. 

If indeed the decoupling condition holds, as discussed, we propose that the problem be 
solved in two stages as illustrated in Fig.l. In the first stage, the MDO approach will be 
applied on a coarse model. The second stage starts with the solution of the MDO algorithm 
and the refined problem is solved with the disciplinary algorithm, thus avoiding the enormous 
complexity of the MDO algorithm when applied on the fine scale model. We claim that the 
resulting design will be a good approximation of the optimal solution. We emphasize that 
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this is possible due to the loose coupling between the two disciplines, otherwise the proposed 
approach will result in a poor approximation of the optimal design since the result will retain 
high residuals of the multidisciplinary optimality conditions. In that case the MDO approach 
should be applied also on fine scales and a preconditioner should be used, as illustrated in 
Eq.(5.5), to achieve effective convergence. 

Finally, how far can we extrapolate the conclusions from this model problem to a more 
realistic model? 

As for the aerodynamic model, it is shown in [12] that an identical symbol for the aerody- 
namic part of the Hessian is obtained when using Euler equations instead of the full potential. 
The analysis for the Navier-Stokes equations has not been completed yet. Shocks were also 
neglected in the aerodynamic model, but we postulate that they are not going to change the 
main conclusion since shocks have a global effect and are not likely to affect the conditioning 
of the Hessian. 

As for the specific modeling which we have chosen to analyze, since there are many 
different models for the cost function and different constraints depending on the application, 
it is impractical to analyze them all. In the model discussed in this paper we used a penalty 
term in the cost function to account for constraints on the structure deformation. However, 
if using inequality constraints instead of penalty terms it is not clear how the coupling of the 
two disciplines will be affected. In practice most of the constraints are not binding at the 
solution, and therefore are effectively of small number. When introduced in small numbers, 
we anticipate that they are not going to change the basic structure of the Hessian near the 
minimum. 
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COARSE MODEL ' H 



Figure 1: An optimization strategy to solve aeroelastic optimization problems in case of 
decoupling as defined in Sec. 5. Apply the MDO approach on a coarse model followed by a 
disciplinary serial approach on fine scales. The result should be a good approximation of the 
multidisciplinary optimal solution. 
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